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ABSTRACT 



We develop an algorithm to reconstruct the triaxial shapes of dark matter 
halos from the anisotropic spatial distributions of their substructures for the con- 
cordance background cosmology. First, we construct an analytic model for the 
anisotropic spatial distribution of dark halo substructures under the assumption 
that the tidal field with non-zero trace in the triaxial mass distribution of the host 
halos generates the substructure bulk motions toward the major principal axes 
of the hosts. Our analytic model implies that the degree of anisotropy depends 
sensitively on the triaxiality of the host halos as well as the correlation between 
the substructure locations and the tidal shear field. Second, we set the axis ratios 
of the triaxial host halos as free parameters in the analytic model, fit the model 
to the numerical results from high-resolution N-body simulation of the concor- 
dance cosmology, and reconstruct the two axis-ratios of the host halos from the 
best-fit values of the free parameters. The comparison of the reconstructed axis 
ratios with the numerical results reveals a good agreement. Finally, we conclude 
that our analytic model may provide a physical understanding of the anisotropic 
spatial distribution of dark halo substructures and a new way to reconstruct in 
principle the triaxial shapes of dark matter halos from the observables. 
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1. INTRODUCTION 

The phenomenon that the satelhte galaxies exhibit strong degree of anisotropy in their 
spatial distributions has recently provoked two-fold controversy: (i) whether the anisotropy 
is along the minor or major axes of the host galaxies; (ii) whether it is consistent with the 
cold dark matter (CDM) paradigm or not. 

In fact, the first observational report of this phenomenon dates back to more than 
three decades ago. Holmberg (1969) analyzed the two dimensional projected positions of 
the satellite galaxies of nearby spirals, and found an effect (the Holmberg effect) that the 
satellite locations within the projected radii < 50 kpc are quite biased toward the regions 
around the minor axes of their hosts. In the beginning, the Holmberg effect was suspected 
not to be real, but due to either the poor number statistics or the projection effects (Hawley 
& Peebles 1975; Sharp et al. 1979; MacGillivray et al. 1982). It is only very recently when 
this effect gains supporting evidence for its existence. Sales & Lambas (2004) used a large 
sample of spiral galaxies and their satellites collected from the Two Degree Field Galaxy 
Redshift Survey (CoUess et al. 2001, 2dFGRS), and found that the satellites have indeed a 
strong tendency to be located near the minor axes of their hosts as Holmberg claimed more 
than three decades ago. 

In spite of the compelling evidence presented in Sales & Lambas (2004), the Holmberg 
effect has yet to be confirmed, because another recent report arouse a new suspicion. Brainerd 
(2005) conducted an analysis of the satellite galaxies about isolated host galaxies in the 
second data release of the Sloan Digital Sky Survey (Strauss et al. 2002, SDSS), and found 
that the anisotropy of the satellite spatial distribution is along not the minor but the major 
axis of the host galaxies, in direct contrast with the Holmberg effect. 

A cosmological turmoil surrounding the satellite anisotropy was caused by the other 
very recent report. Kroupa, Theis, & Boily (2005) showed that the Milky Way (MW) dwarf 
satellites are located near polar plane (Lynden-Bell 1982; Majewki 1994), and claimed that 
this anisotropic distribution of MW satellites is inconsistent with being drawn from a CDM 
substructiu^c population. Shortly after Kroupa et al. (2005), a fiurry of numerical studies 
demonstrated that the CDM model indeed predicts anisotropic satellite distributions. Kang 
et al. (2005) showed by their high-resolution N-body simulation that the satellites are likely 
to lie in the plane defined by the major and intermediate axes of the host halo. Zentner et al. 
(2005) confirmed the findings of Kang et al. (2005) by providing the numerical results that 
the dark halo substructures on galactic scales are preferentially located along the major axes 
of the triaxial mass distributions of their hosts, proving that the CDM paradigm does not 
predict the isotropic spatial distribution of sateUite galaxies (see also, Knebe et al. 2004). 
Libeskind et al. (2005) also performed high-resolution simulations of the galactic CDM halos 
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to verify that the sateUite anisotropy is not inconsistent with the CDM cosmology. 

Yet, an important problem highlighted by Kroupa et al. (2005) still remains a mystery: 
why are those halos forming the luminous dSphs located in a disk-like configuration almost 
perpendicular to its host galaxy ? This question in fact holds a crucial key to understanding 
the formation of galaxies since if dSphs are assumed to follow the dark matter distribution, 
the observed polar alignment of dSphs indicates a disagreement between the baryons of the 
host galaxy and its dark matter halos. 

To address this issue and confirm the Holmberg effect, it is necessary first to under- 
stand the physical origin of the satellite anisotropy. Now that several numerical studies 
proved that the satellite anisotropy is consistent with the CDM model, it should be pos- 
sible to construct a theoretical model for the anisotropic spatial distribution of satellites 
in the CDM context. The standard scenario explains that the anisotropy in the satellite 
distribution is originated from the anisotropic infall of material along filaments (West 1989, 
Bond, Kofman, & Pogosyan 1996). But, their explanations were inconclusive, making no 
quantitative prediction for the degree of anisotropy expected on different scales. 

Very recently, Lee, Kang, & Jing (2005, hereafter LKJ05) pointed out the incompleteness 
of the anisotropic infall scenario in their studies of the anisotropic orientations (not the 
anisotropic locations) of dark halo substructures. They claimed that since the anisotropic 
infall is a primordial effect, being quickly damped away by the subsequent nonlinear process 
like the secondary infall, it is unlikely to be the direct cause for the anisotropic orientations of 
substructures. They assumed that the substructure anisotropic orientations must be caused 
by the tidal interactions with the gravitational fields in the triaxial mass distribution of their 
host halos, and constructed an analytic model for the alignments between the major axes of 
the substructures and their host halos from the physical principles, which turned out to be 
in good agreement with simulation results. 

Motivated by the success of LKJ05 model, here we attempt to construct a similar 
analytic model for the anisotropic spatial distribution of dark halo substructures assuming 
that the anisotropy in the spatial locations of substructures are also induced by the tidal 
interaction with the gravitational field in the triaxial mass distribution of host halos. Our 
final goal is to develop an algorithm for reconstructing the triaxial shapes of host halos from 
the observable anisotropic spatial distribution of their substructures by using this analytic 
model. Throughout this paper, we assume a concordance ACDM cosmology with = 0.3, 
^A,o = 0.7, and h = 0.7. 
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2. PHYSICAL ANALYSIS 

Lee & Pen (2001) claimed that the position vectors of dark halos are related to the tidal 
shear tensors, and suggested that the relationship between the two quantities can be written 
as a first-order correlation (see eq. [12] in Lee & Pen 2001). They supported their claims 
by N-body simulations (see also Lee & Pen 2000) where it was found that the tidal shear 
tensors are indeed correlated with the halo position vectors at first order. Although they did 
not show that it also holds on the substructure scale, we extrapolate here their claim to the 
substructure scale, and adopt an Ansatz that the position vectors of dark halo substructures 
are correlated with the tidal shear tensors of their host halos at first order. 

Suppose that a subhalo lies in the perturbation potential field inside its triaxial host 
halo. The potential ^ is assumed to be smoothed on the host halo mass scale. Let q = (qi) 
be the initial position of the subhalo when there is no interaction with the host halo potential 
field, and let also x = (xj) represent the final position of the subhalo under the interaction 
of the host halo potential. Both vectors, x and q, are defined in the center of mass frame of 
the host halo. Following the linear scheme (Zeldovich 1970), the relation between (xj) and 
(qi) in the host halo potential may be written as 

Xi = -D+di-^{q), (1) 

where is the linear growth rate of the density field. Note that x in equation (1) represents 
the displacement of the subhalo position generated by the interaction with their host halo 
gravitational field. 

Taylor-expanding ^ to second order around the host halo center of mass, we have 

^ = *(0) + qkdj^lo + ^qjqkdjdk^lo. (2) 
Differentiating equation (2) gives 

di-^ ^di-^\o + qkdidk-^\o. (3) 
Putting equation (3) into equation (1), we end up with 

Xi oc qkTik{0), (4) 

where T = (Tj^) is the tidal shear tensor defined as didk^ . Equation (4) implies that the 
subhalo displacement is a bulk motion generated by the nonzero quadrupole moment of the 
triaxial gravitational potential of the host halo. 

One might think that since equation (1) is derived from the linear perturbation theory 
based on the Zel'dovich approximation, it should be invalid in the highly nonlinear regime. 
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But, Lee & Pen (2001) showed by N-body simulations that even in the highly nonlinear 
regime the correlation between the halo position vector and the tidal shear tensor can be 
written as first order like equation (4) 

To investigate how the direction of x is distributed given T, we first consider the co- 
variance {xiXj\T). Using equation (4), we have {xiXj\T) oc (qkqiTikTij) . If the two tensors, 
(qiqj) and {TikT^j), were mutually independent, then, wc would end up with having a simple 
expression, {quqiTikT^,) = {qkqi)TikTij oc TikT^j as {quqi) oc 5ki. However, it it not necessarily 
true that {qiqj) is uncorrelated with {Tn-T^j). 

To take into consideration the possible correlation between (qiqj) and (TikTkj), and to 
overcome the limited validity of the linear approximation, we use the following parametrized 
formula as Lee & Pen (2001) did: 

I — s ^ ^ 

{xiXj I T) = -^^ij + sTikTkj , (5) 

where Xi is the rcscalcd position vector (rcscalcd in a sense that it docs not have a dimension 
of "length"), T = (T^j) is the unit tidal shear tensor with non-zero trace: Tik = ^iA;/|T| 
and Tr(T) ^ (we use the notation of T to distinguish it from the traceless unit tidal 
shear tensor T), and s is a free parameter in the range of [—1, 1] to quantify the correlation 
between (xiXj) and {TikTkj)- In the linear theory, the value of s is given as unity. 

Note that equation (5), although it is derived from the linear perturbation theory, has 
already included the nonlinear effect. The correlation parameter, s, quantifies how the 
correlation between x and T differs from the linear theory prediction. The degree of the 
deviation of its value from unity measures how strong the nonlinear effect is. 

It is also worth noting that equation (5) bears a close resemblance to the quadratic 
formula used by LKJ05 for the correlation between the dark halo angular momentum L = 
(Lj) and the unit traceless tidal shear field T = (Tj^) defined as T^j = Tij/\T\ with = 
%j - 'IV(T)5,,/3 

1 + c - - 

(LjLj|T) = ~^^i3 ~ cTikTkj, (6) 

where the correlation parameter c is in the range of [0,1]. LKJ05 used equation (6) to 
investigate the alignments between the major axes of substructures and their host halos. 

Similar as equations (5) and (6) look, it is very important to recognize the key differences 
between the two. In equation (6), T is traceless. Consequently, L is preferentially aligned 
only with the intermediate principal axis of T. Whereas in equation (5) T has non-zero 
trace. Consequently, x is preferentially aligned with either the minor or the major principal 
axis of T: If s is negative, then x is aligned with the minor axis of T. If s is positive, x is 
aligned with the major axis of T. 
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Strictly speaking, equation (5) holds only if T and x are smoothed on the same scale. 
Here, they don't: T is supposed to be smoothed on the host halo mass scale, while x on the 
substructure radii. For simplicity, here, we just assume that equation (5) still holds, ignoring 
the difference between the host halo and their substructure smoothing scales. 

In a similar manner with LKJ05, assuming that the probability density distribution of 
p(x|T) is Gaussian, and expressing the position vectors x in terms of the spherical polar 
coordinates in the principal axis frame of T as x = (a; sin a cos x sin a sin /5, x cos a) where 
X = |x|, q; and (5: the polar and the azimuthal angles of x, respectively, we derive the 
probability density distribution of the cosines of the polar angles, p(cosQ;) as: 



p(cos a) = ^ n (l - 5 + 3sA ■ ) ^ X 



=1 

3 

sin^ a cos^ (5 sin^ ol sin^ /9 cos^ a ^ ^ 



Vl-s + 3sA? l-s + 3sA^ l-s + 3sAi 



(7) 



A? + A^ + A^ = l, Ai + A, + A3= + + (8) 

where {A^jf^Q arc the three eigenvalues of the tidal shear tensor T in a decreasing order, and 
{Ai}f=i are the unit eigenvalues of T defined as Aj = Aj/ (A^ + A| + A|)^^^. 

In the hnear density field, the sum of all the three eigenvalues of the tidal shear tensor 
equals the linear density contrast: Ai + A2 + A3 = 5. According to the spherical dynamical 
collapse model, a gravitational bound object forms when the linear density contrast reaches 
its critical value, be- Hence, one can say that inside the host halo Ai + A2 + A3 = b^. The 
value of the critical density contrast, ^c, scales with the redshift z as bjyz) = Dj^{z)6qc where 
6co is the critical value at 2; = 0, which depends very weakly on cosmology (Eke, Cole, & 
Frenk 1996). For example, Sco = 1.686 for the flat universe of closure density. Now, one can 
rewrite equation (8) as 

Ai + A2 + A3 = 77^ (9) 

(A? + Ai + Ai)^/^ 

As LKJ05 did, wc assume further that the tidal shear and the inertia shape tensors of the 
host halos are strongly correlated with each other in an opposite direction (Lee & Pen 2000; 
Porciani, Dekel, & Hoffman 2002), so that the minor and the major axes of the tidal shear 
tensor are the major and minor axes of the inertia shape tensor. Under this assumption, 
equation (7) basically represents the probability density distribution of the cosines of the 
angles between the position vectors of the substructures and the major principal axes of 
their host halo. Equation (7) reveals that p(cosQ;) depends sensitively on the value of s: If 
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— l<s<0,xis aligned with the major axis of the host halo (the minor axis of T); If 
< s < 1, X is aligned with the minor axis of the host halo (the major axis of the T); If 
s — 0, there is no alignment between x and any principal axis of the host halo. 

Note that even when s 7^ 0, it is possible to have no alignment between x and the host 
halo axes. The distribution p(cos«) depends not only on the value of s but also on the unit 

eigenvalues with non-zero trace, {Ajjf^Q. For the case of the traceless eigenvalues of {\i}f^Q 
which satisfies the two constraints: Yl^=o^i ~ ^ ^^'^ SLo ~ ^ (^^^ §2 LKJ05), their 
values are well approximated as the average Ai = —A3 = 1/V2 and A2 = regardless of 
the original values of {\i}f^i. However, for {Ajjf^Q, the constraints (cq. [8]) arc expressed in 
terms of the original values of {Ajj^^Q. Therefore, even when the correlation parameter s 
has the maximum absolute value (i.e., |s| = 1), if the host halo has almost spherical shape 
satisfying Ai = A2 = A3 = 0, then p(cos a) becomes uniform, implying no anisotropy in the 
spatial distribution of substructures. 



3. HALO SHAPE RECONSTRUCTION 

To find p{cos a) numerically, we use the same numerical data as used in LKJ05 which 
were generated by P^M N-body simulations with 256^ particles in a 100/i~^Mpc cube (Jing 
& Suto 1998, 2000) for the ACDM cosmology with = 0.3, Qa,o = 0.7, and h = 0.7. In 
the simulation data, we selected four dark halos on mass scales around 5 — 10 x 10^^/i~^Mq, 
and identified the self-bound substructures in each host halo by using the SUBFIND routine 
of Springel et al. (2001). For the detailed descriptions of the N-body simulations and the 
subhalo identification procedure, see §3 in LKJ05. 

For each host halo selected in the numerical data, we compute the inertia shape tensor 
lij = ^a''^afa,ifa,3 ^ud determine the major principal axis by finding the eigenvector eh 
of (Jjj) corresponding to the largest eigenvalue. In each host halo principal axis frame, 
we measure the unit position vectors of the subhalos x^, and the cosines of the angles, a, 
between the major axis of the host halo and the locations of their subhalos by computing 
cos a = |e/i ■ Xs|. We find the probability density distribution, p(cosQ;) by counting the 
subhalo number density as a function of cos a. 

We fit our analytic model (eq.[7]) to this numerically determined p{cosa) with setting 
the three parameters, Ai,A2 (A3 is not independent according to eq.[8]) and s as free, and 
determine the best-fit values of the free parameters by means of the statistics. We perform 
this procedure for each individual host halo (CL#1, CL#2, CL#3, CL#4) separately at four 
different redshifts {z — 0, 0.5, 1, 1.5). 
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Figure 1 plots the final fitting results. In each panel, the histogram represents the 
numerically measured p(cosQ;), while the (red) solid curve corresponds to the analytic dis- 
tribution (eq. [7]) with the best-fit correlation parameter, s. As one can see, for all cases, 
the best-fit correlation parameter s has a negative value. In other words, the dark halo 
substructures are preferentially located near the major axes of their host halos, consistent 
with previous numerical findings (Knebe et al. 2004; Zentner et al. 2005; Libeskind et al. 
2005). Indeed, the analytic distributions agree with the numerical results for all four clusters 
at all redshifts quite well. 

As a final step, wc reconstruct the triaxial shapes of the four host halos from their 
distributions p(cosQ;). For this purpose, we adopt the following formula suggested recently 
by Lee, Jing & Suto (2005): 



02 /l-D+AaV^' «3 /1-I?+Ai^'/' 



(10) 



where oi, 02, 03 represent the three axis lengths of the triaxial host halos in a decreasing order. 
Equation (10) basically links the two axes ratios of the triaxial host halos to the eigenvalues 
of the tidal shear tensors. According to this equation, once {Aj}|^3 is determined, then one 
can reconstruct the axis-ratios of the host halos. 

Equations (8) and (9) allow one to determine the values of {Aj}f^3 from the best-fit 
values of {Ai}f^3 as 

, (^cAi 4A2 . 4 A3 , , 

^ ~ 7 ? ^^2 — Y r ^3 — ^ [i-i-) 

Al -|- A2 -|- A3 Al + A2 + A3 Al + A2 + A3 

We use Sc = 1.686 as an approximation and the best-fit values of Ai, A2, A3 to determine the 
values of Ai,A2,A3 using equation (11) for each host halo, and finally reconstruct the axis- 
ratios of the host halo using equation (10). With the above method, we finally reconstruct 
the two axis-ratios of all four clusters at four redshifts. 

Figure 2 plots the intcrmcdiatc-to-major axis ratio vs. the minor-to-major axis ratios 
of the four clusters at four redshifts. In each panel, the (red) triangle dot represents the 
reconstructed axis-ratios by the above algorithm, and the (blue) square dot is the original 
axis-ratios from simulations. Although the reconstruction does not look perfect, it works 
quite satisfactorily, especially given all the approximations used in the algorithm. It is 
obvious that at least one of the two axis-ratios is fairly accurately reconstructed. We note 
that for the case of s ~ —1, the difference between 02/01 and O2/01 is relatively large, which 
trend is also reconstructed well by our algorithm. 
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4. DISCUSSIONS AND CONCLUSIONS 

Our achievements are summarized as follows: (i) we construct an analytic model for 
the anisotropic spatial distribution of dark halo substructures; (ii) we confirm numerically 
that the dark halo substructures have a strong propensity to be located near the major axes 
of their host halos; (iii) we develop an algorithm to reconstruct the triaxial shapes of dark 
halos from the anisotropic spatial distribution of their substructures. 

Successful as are our analytic model and the reconstruction algorithm based on it, it is 
worth noting a couple of its inconsistencies. First, as mentioned in §2, equation (7) holds only 
when the tidal shear field is smoothed on the substructure radii. We ignored the smoothing 
scale difference, which should cause some scatters in the measurement of the correlation 
parameter s. Second, we adopted the spherical collapse condition Ai + A2 + A3 = 5c ~ 1.686 
even though we consider triaxial dark matter halos. Definitely, it will be necessary to refine 
our analytic model further, accounting for these two inconsistencies. 

Another limitation of our model is that its underlying assumption has been numerically 
tested only for the case of the concordance cosmology. Equation (10) was shown to work only 
provided that the background cosmology is described by the concordance cosmology (Lee, 
Jing & Suto 2005). To find the dependence of the correlation parameter, s, on cosmology, 
and to apply the model to a broader range of cosmological models, it will be first necessary 
to make more general assumptions and to test the underlying assumptions against numerical 
results, which is definitely beyond the scope of this paper. 

Now that we proved that it can reconstruct in principle the axis-ratios of triaxial dark 
halos from the anisotropic spatial distribution of substructures for the concordance ACDM 
model, it is worth discussing about how to apply the algorithm to real data from obser- 
vations of such as the weak gravitational lensing. In our analytic model, the substructure 
positions are defined in the real three dimensional space relative to the principal axes of host 
halos, which is, however, very difficult to determine in practice. In observations they are 
measured in the two dimensional space projected along the line-of-sight direction which is 
not necessarily coincident with the host halo principal axis. Therefore, to apply our recon- 
struction algorithm to real data, it will be necessary first to find an analytic expression of 
the substructure spatial distributions in the general non-principal axis frame and to account 
for the effect of projection along the line-of-sight. Our future work is in this direction, and 
we wish to report it elsewhere soon. 

Another direction is to explain the Holmberg effect. As shown, our analytic model 
for the substructure spatial distribution (eq. [7]) is characterized by the single correlation 
parameter, s. By investigating what physical mechanism caused the correlation parameter 
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to change their signs, it might be possible to explain the Holmberg effect theoretically. 

Finally, we conclude that our analytic model may provide a physical quantitative un- 
derstanding of the anisotropic spatial distribution of dark halo substructures, and will allow 
us to determine in principle the axis-ratios of triaxial dark halos from the observables. 
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Foundation. X. K. is partly supported by NKBRSF (G19990754), by NSFC (Nos. 10125314, 
10373012), and by Shanghai Key Projects in Basic Research (No. 04jcl4079) 
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Fig. 1. — Probability density distributions of the angles between the major axes of four 
individual host halos and the position vectors of their substructures at four different epochs; 
z = 0, 0.5, 1 and 1.5. In each panel, the histogram represents the simulation results, and the 
(red) solid line represent the theoretical prediction (7). 
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Fig. 2. — Axis-ratios of four individual clusters at four different epochs; In each panel, the 
square represents the real axis-ratio measured in simulations while the triangle represents 
the reconstructed axis-ratio. 



